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Abstract 



Using molecular dynamics calculations and the Voronoi tessellation, we 
study the evolution of the local structure of a soft-sphere glass versus temper- 
ature starting from the liquid phase at different quenching rates. This study 
is done for different sizes and for two different boundary conditions namely 
the usual cubic periodic boundary conditions and the isotropic hyperspherical 
boundary conditions for which the particles evolve on the surface of a hyper- 
sphere in four dimensions. 

Our results show that for small system sizes, crystallization can indeed be 
induced by the cubic boundary conditions. On the other hand we show that 
finite size effects are more pronounced on the hypersphere and that crystal- 
lization is artificially inhibited even for large system sizes. 
PACS numbers: 61.43.Bn,64.70.Pf,61.50.Ah 
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When cooling a liquid down to low temperatures, what is the nature of the resulting solid 
phase ? The answer to this question in experimental studies or in numerical studies depends 
on how fast the liquid has been cooled down i.e. on the quenching rate. If the liquid is cooled 
faster than a critical rate, the resulting phase will be a glass otherwise the system will evolve 
towards the more stable crystal phase. The knowledge of this critical rate is very important 
especially in numerical studies of model glasses since crystallization effects should not bias 
the results at low temperature. Similarly to what had been done in a Lennard- Jones system 
IjH or in liquid sodium |2| we used the Voronoi cell statistics to determine this critical rate 
in a model soft-sphere glass . Nevertheless in numerical studies due to computer time and 
memory limitations, one has to mimic macroscopic samples with the use of the so-called 
"boundary conditions". This is usually done by replicating a cubic box containing a few 
thousands of particles, for which the calculations will be explicitly done, in the 3 directions 
of space 0. These boundary conditions are generally called periodic boundary conditions 
(PBC) and let us call this periodic space C3. Obviously in this geometry the finite systems 
stay homogeneous but can bear a slight anisotropy. Therefore it is justified to address the 
question of the influence of this anisotropy on the crystallization of our model glass for which 
the ground state is a face-centered-cubic arrangement. In other words, do the cubic PBC 
favor crystallization ? 

In this letter we investigate the influence of the boundary conditions and the system size 
on the evolution of the structure of a soft-sphere system during the quench with the use 
of classical molecular dynamics (MD) simulations. To distinguish an amorphous structure 
from a crystal or a partially crystallized system, we use the fraction of pentagonal faces of 
the Voronoi cell attached to each particle, which is an excellent tool to obtain this distinction 
0. As an alternative to the cubic PBC we use the "hypersphere boundary conditions" . 
This means that the MD simulations were performed by confining the system of soft-spheres 
to the surface 5*3 of a hypersphere in four dimensions. This technique was used earlier for 
systems of charged particles |5| , or fluids of hard-spheres || and Lennard- Jones particles [|7| . 
The space 53 (as well as C3) is homogeneous, finite and non-Euclidean; but in addition, S3 
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is isotropic therefore there is no preferential direction which could influence crystallization. 
Moreover for the calculation of long-range forces like Coulomb forces an advantage of using 
5*3 is a reduction of computer time by a factor 3-4 compared to the standard Ewald method 
||. Of course asymptotically for very large system sizes the results in C3 and 5*3 should 
converge. 

Our results show that for small system sizes crystal nucleation occurs in C3 while the system 
remains amorphous in S3. This is coherent with previous results obtained for hard-sphere 
collections |J and shows that the glass stability is, as desired, higher on the hypersphere. 
Nevertheless for larger systems we show that, compared to the asymptotic behavior which 
is as expected identical in C3 and S3 , the use of hyperspherical boundary conditions can 
artificially discourage crystal nucleation and therefore the value of the critical quenching 
rate can be underestimated. These results permit to shed new light on the influence of the 
widely used PBC on crystallization effects in model glasses and show the consequences of 
using isotropic boundary conditions in a curved three dimensional space such as 5*3. 

In our microcanonical MD simulations the particles interact via the repulsive potential 
proposed by Laird and Schober ||: 

U(r) = e(^) 6 + Ar 4 + B for r < 3a (1) 
= for r > 3cr 

The constants A and B are chosen so that both the potential and the force are zero at the 
cut-off distance 3a. To avoid reduced units, we use the values of Lennard- Jones argon for 



the mass m of the particles, as well as for e and a |TU| . To integrate the equations of motion 
we use the Velocity Verlet algorithm ||] with a timestep At = 2.5 fs. To investigate the 
influence of the system size we use samples containing N = 512, 1000 and 5000 particles and 
to get better statistics we use either 10 independent configurations for iV = 512 and 1000 or 2 
configurations for N = 5000. For each system we start from a well equilibrated liquid sample 



at a temperature around 50K (well above the melting temperature T m m 23K |TTJ) which is 



then cooled down to zero temperature at two different quenching rates: 4.0 x 10 11 K/s and 
10 11 K/s. The fastest quenching rate is close to the critical rate for which no crystallization 
should be observed after the quench while the second rate is smaller then this critical rate and 
therefore the system should be in a crystallized or partially crystallized state at the end of 
the cooling process || . To study the evolution of the structure as a function of temperature 
we save all along the quenching procedure configurations (positions and velocities) on which 
the Voronoi tessellation scheme is applied (it is worth noticing that no relaxation period has 
been considered). Once the Voronoi cell attached to each particle has been determined, we 
can calculate the fraction f$ of pentagonal faces in the system. This is a good indicator of 
the degree of randomness (or on the other hand the degree of crystallization) since a large 
value of fs (typically around 0.45) is a sign of strong icosahedral local order characteristic 
of a glass phase, while a small value of /s (< 0.2) indicates crystal nucleation. Therefore we 
know for each sample if crystallization took place during the quench and we can determine 
at which temperature this crystallization started. 

While MD simulations in C3 are common and well established, we think it useful to give 
more details on the simulations in S3. First of all concerning the density p, it is fixed and 
we use the value p = Na 3 /V space = 1. In C 3 we use standard PBC with a rigid cubic box 
of edge length L (y space = L 3 ) while in S3 V space = 2n 2 R 3 where R is the radius of the 
hypersphere. Fundamentally the volume fractions of the corresponding hard-sphere systems 
are not strictly equal in both spaces because one should take into account the curvature of 
the hard-spheres in the fourth dimension. Nevertheless a straightforward calculation shows 
that the relative error in S3 is smaller than 5 x 10~ 3 % for the smallest system considered 
here (N = 512). Second let us detail a MD step in S3. For convenience we chose to 
use the quaternion formalism [HJ] in order to describe the equations of motion in the four 
dimensional space. The unit quaternion Q l = (q^, q\, q\, q\) attached to particle i is defined 
by the set of q l k s given by q\ = x\/R, where the x k s are the coordinates of particle i in four 
dimensions, verifying Y2k ( x k) 2 = R 2 - We consider the potential U(r) as being a function of 
the interparticle cartesian distance r in the four dimensional space. Thus the force F l on 
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particle % is a four dimensional vector, defined by the components f k given by: 

i no- 1 7v W 



where \Q^\ is related to the distance between particles i and j by = r^jR. Next 
the force is projected on the three dimensional tangent space of particle i by calculating the 



quaternion product 12 



¥ = Q i F i (3) 

where Q l = (q^, —q\, —q l 2 , —q\) is the "conjugate" of Q\ The acceleration of particle i in the 
three dimensional tangent space Y k is obtained using the Newton's equations, 

f k = n/m , k = 1,2,3 (4) 

Once the acceleration has been calculated, the velocity is recovered using the standard 
Velocity Verlet scheme ||] and then the particle i can be moved to its new positions Q 1 ' 
using 

n* 6QlQl ™ 

with 5Q i = (l,6q\,6ql,6q$) and 5g£ = (Atvj. + (At 2 /2)f k )/R. The determination of the 
new positions completes the iteration and the process can be repeated as desired. In the 
scheme described above the dynamics is in fact performed in the four dimensional space with 
the constrain that the points should stay on the hypersphere. This implies that the first 
component of the force /q, orthogonal to the hypersphere, is compensated by an adequate 
reaction force. We could have considered an alternative procedure consisting in staying on 
the hypersphere and considering forces between particles separated by the "geodesic" dis- 
tance ||. But, since only distances smaller than 3a are considered here (for larger distances 
the potential is zero), the relative difference (due to the curvature of the hypersphere) be- 
tween the cartesian distance and the geodesic distance remains smaller than a fraction of a 
percent for the smallest system considered here (512 particles). Therefore the quantitative 



differences between the two approaches stay very small and in fact vanish in the infinitely 
large system limit. 

Finally concerning the Voronoi tessellation in C3, the scheme has been detailed elsewhere 
|T3|j . In S3 the tessellation is done in the four dimensional space using the geodesic distance 
between particles |T4|]. For all the different samples and quenching rates, all the calculations 
have been made in parallel in S3 and in C3 in order to compare the effect of the boundary 
conditions on the final structure of the samples after the quench. 



First we consider the fastest quenching rate (4.0 x 10 11 K/s). In Fig. 1 we report on the 
variation in C3 and £3 of f$ versus the temperature for the three sizes considered in our 
study (N = 512, 1000, 5000). It is worth recalling that the results have been averaged over 
several samples for the different values of N. The first observation concerns the asymptotic 
behavior obtained for the largest systems in both spaces: as expected this behavior is sim- 
ilar. Indeed if L or R increases, the description becomes closer to the one of the infinite 
system. As a general result, the finite size effects are more important for temperatures below 
the melting temperature. Moreover, if one excludes the case N = 512 in C3, it appears that 
these finite size effects are more important in S3 than in C3. This is clearly visible when 
comparing the curves for N = 5000 and N = 1000 which are almost superimposed in C3 and 
not in 5*3. These finite size corrections act in opposite directions when comparing C3 and 
S3: in C3 the small system sizes favor a small value of f$ while in S3 high values are favored. 
Indeed, the smaller the hypersphere, the higher f 5 i.e. the more the system is icosahedral. 
This can be explained by the fact that for N = 120 the ground state in S3 is the "Polytope 
{3,3,5}" |L5[] made of dodecahedral Voronoi cells only, which means that for this ideal struc- 
ture / 5 = 1. Therefore the more N decreases and comes closer to 120 the more the structure 
will resemble the Polytope {3,3,5} and the more fa will increase and tend towards 1. The 
fact that even for large systems f§ remains greater than 0.43 at T=1K, indicates that no 
crystallization occurred during the quench and shows that the quenching rate is above the 
critical rate separating glass and crystal forming rates. Another point concerns the glass 
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transition which occurs around T g « 10. 5K || and which is reflected in a saturation of f$ 
below T g . This saturation is not obvious for N = 512 but is already visible for N = 1000 
and is much more pronounced for N = 5000. Our results do not show a significant shift of 
T g as a function of system size. 

The huge difference between C3 and S3 is observed for iV = 512. Indeed in C3, for iV = 512, 
Fig. 1 shows clearly a drop of fs below T g with a saturation towards 0.36 at low temper- 
ature. Since these results have been averaged over ten samples, this indicates that some 
of the samples have crystallized during the quench. Snapshots of the structures at T=1K, 
confirm that three out of ten samples are crystals, one of them having a fraction of pentago- 
nal faces smaller than 10~ 2 indicating an almost perfect crystalline character. These results 
are coherent with what has been observed in hard-sphere collections || and indicate that 
for small systems (a few hundred particles) the cubic PBC can indeed induce crystallization 
while the same systems remain "perfectly" amorphous with hyperspherical boundary con- 
ditions. Nevertheless this "encouragement" to crystal nucleation in C3 already disappears 
for systems containing 1000 particles which is the size we have used in earlier studies [Q. 
In a second step we consider a quenching rate of 10 11 K/s. This rate is smaller than the 
critical rate and therefore the structure after the quench should be crystalline. In Fig. 2 we 
show the variation of f$ versus the temperature in 5*3 and C3 for iV = 1000 and iV = 5000. 
As expected the asymptotic behavior (N = 5000) is the same in both spaces: a drop of 
f 5 is observed below T g indicating crystallization. In C 3 this effect is even more spectacu- 
lar for the systems containing a thousand particles. In fact almost all of the ten samples 
considered in this study show a value of /s smaller than 0.2 at T=1K indicating that these 
samples have undergone crystallization. A look at snapshots of these crystalline structures 
does not show an obvious correlation between the axes of the periodic cubic box and the 
crystallographic directions. The fact that the samples containing N = 1000 particles exhibit 
an almost perfect crystalline character may be due to the difficulty to build extended defects 
in such small crystalline systems. This is not the case for the 5000 particles samples which 
may explain why these samples contain more defects and therefore show a higher value of f$ 

7 



at small temperature (/ 5 0.35). Nevertheless crystal nucleation occurred in S 3 and C 3 for 
N = 5000 and in C 3 for N = 1000. On the contrary the systems containing a thousand par- 
ticles in S3 do not show any sign of crystal nucleation. This indicates that for small system 
sizes the hyperspherical boundary conditions discourage artificially crystal nucleation. If 
the aim is to avoid crystallization then one should indeed work in S3 with a relatively small 
system. On the other hand if one wants to determine the critical rate separating crystal 
forming rates from glass forming rates then the use of hyperspherical boundary conditions 
can lead to an underestimation of this rate. Indeed we quenched liquid samples containing 
1000 particles at a rate of 4 x 10 10 K/s and still we did not detect any sign of crystallization 
in S 3 even though the 5000-particles system crystallized. This is another example of the 
strong size effects in S 3 : one needs to use a large system to get the correct results. On the 
contrary in C 3 the correct behavior is already obtained for systems containing 1000 particles. 

In conclusion, we have used molecular dynamics simulations and the Voronoi tessella- 
tion to study the influence of the boundary conditions on the structure of the solid phase 
obtained after the quench of a soft-sphere system. This has been done for several system 
sizes (512, 1000 and 5000 particles) and for two different quenching rates close to the critical 
rate separating crystal from glass forming rates. To simulate a macroscopic sample we have 
used the usual cubic periodic boundary conditions (PBC) but also isotropic hyperspherical 
boundary conditions for which the particles are confined to the surface of a hypersphere in 
four dimensions. For the glass forming rate we show that for the smallest system size, the 
PBC can indeed induce crystal nucleation while the larger samples remain amorphous. This 
shows that caution should be used with the usual PBC when samples containing only a few 
hundred particles are considered. On the contrary for the crystal forming rate the results 
show that the hyperspherical boundary conditions can discourage artificially crystal nucle- 
ation even for samples containing thousand particles. On the one hand this is an advantage 
because no crystallization effects bias the results supposed to be obtained in a glass. On the 
other hand for a given system it leads to an underestimation of the critical quenching rate. 
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Of course the asymptotic results obtained with 5000 particles are as expected similar using 
both boundary conditions but since it appears that the finite size effects are smaller when 
the usual PBC are used and even though the calculations on the hypersphere are slightly 
faster, our results indicate that in any case the correct behavior is obtained with the usual 
PBC for the 1000-particles samples. Of course this critical size depends on the interaction 
potential and should be re-determined for each different system. 
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FIGURES 

FIG. 1. Fraction of pentagonal faces, /s, versus temperature in C3 (open symbols) and 53 
(black symbols) for N = 512 (triangles), N = 1000 (squares) and N = 5000 (circles). 
The quenching rate is equal to 4 x 10 11 K/s. 

FIG. 2. Fraction of pentagonal faces, fa, versus temperature in C3: N = 1000 (A) 
and N = 5000 (•). Idem in S 3 : N = 1000 (♦) and N = 5000 (■). 
The quenching rate is equal to 10 11 K/s. 
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Fig. 2 
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